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Abstract 

We show that the stochastic Morris-Lecar neuron, in a neighborhood of its stable point, 
can be approximated by a two-dimensional Omstein-Uhlenbeck (OU) modulation of 
a constant circular motion. The associated radial OU process is an example of a leaky 
integrate-and-fire (LIE) model prior to firing. A new model constructed from a radial OU 
process together with a simple firing mechanism based on detailed Morris-Lecar firing 
statistics reproduces the Morris-Lecar Interspike Interval (ISI) distribution, and has the 
computational advantages of a LIE. The result justifies the large amount of attention paid 
to the LIE models. 
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1. Introduction 

Much effort has been made to create a realistic but still easily computed stochastic neuron 
model, primarily by combining subthreshold dynamics with firing rules. The result has been 
a variety of, usually one dimensional, leaky integrate-and-fire (LIF) descriptions with a fixed 
membrane potential firing threshold H [TT] [TSl [19], or with a rate of firing depending more 
sensitively on membrane potential ifTSl |2TI . These models are useful both for obtaining 
analytical results and for ease of simulation. 

By contrast, the two-dimensional stochastic Morris-Lecar (ML) neuron model, a simple 
cousin to the more detailed Hodgkin-Huxley (HH) model, describes the dynamics of firing in 
a way more closely motivated by the biology. It has been better respected by biologists than 
the LIF class of models, but has received little attention owing to the difficulty of mathematical 
analysis of this rather complicated stochastic dynamical system. 

In Section|4]of this paper we show that in fact a LIF model is embedded in the ML model as 
an integral part of it, closely approximating the subthreshold fluctuations of the ML dynamics. 
This result suggests that perhaps the firing pattern of a stochastic ML can be recreated using 
the embedded LIF together with a ML stochastic firing mechanism. We construct such a model 
in Section|5]and|6j and show in Section|7]that its Interspike Interval (ISI) distribution is similar 
to that of the ML. Our model, while of the type described in our first paragraph, combines the 
realism of the ML with the ease of analysis and computation of a one dimensional LIF-type 
model. The work invested in LIF models is further justified by this new model. 
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Before we set up our stochastic ML model and write analytical details, let us have an 
informal look at how it works. The principal dynamics of the ML, in the central range of the 
input current, consist of a stable limit cycle (Fig. [T]^) corresponding to firing, which encloses 
a stable fixed point. In between there loops an unstable limit cycle. The path of the stochastic 
model has two quasi-stable patterns (Fig. [TJj). One is succesive firings, where the dynamics 
makes "large" noisy circuits around the stable limit cycle, the other is membrane fluctuations 
between spikes, where the dynamics makes "small" noisy circuits around the fixed point inside 
the unstable limit cycle. The system would continue forever in one of these two patterns were 
it not for the noise which causes switching from firing to subthreshold fluctuations and back 
again at random times when the dynamics cross the unstable limit cycle. Our analysis will 
show that the dynamics between spikes, of random cycling inside the unstable limit cycle 
followed by crossing to the stable limit cycle outside it, can be identified with the sample path 
behavior of a two-dimensional Ornstein-Uhlenbeck (OU) process times a rotation. 

A main ingredient in our result is the stochastic dynamical phenomenon that oscillations 
which damp to a fixed point in a deterministic system will be sustained by the stochasticity 
in a corresponding stochastic system. Damped oscillations in a two-dimensional system are 
signalled by a local linear structure defined by a matrix having a pair of conjugate complex 
eigenvalues with negative real part. A corresponding stochastic system will not damp, being 
prevented by the noise. Instead, a quasi-stationary stochastic process is set up, which cycles 
in a random pattern around the fixed point. Using recent results of [2 1 we are able to identify, 
approximately, this stochastic process which is part of the subthreshold dynamics of the ML. 
Up to a fixed linear transformation, the approximating process is the product of a steady fast 
rotation with a two-dimensional OU process. The identification allows us to cement in place 
the correspondance, for a particular set of model parameters, a particular LIF model as the 
appropriate subthreshold phase between ML firings. 

2. The Morris-Lecar model 

There exists a large variety of modeling approaches to the generation of spike trains in 
neurons (see e.g. ||6l[TT][T4l). Most famous is the Hodgkin-Huxley (HH) model [ 13 1 consisting 
of four coupled differential equations, one for the membrane voltage, and three equations 
describing the gating variables that model the voltage-dependent sodium and potassium chan- 
nels. A large amount of research effort is currently directed towards understanding how neural 
coding carries information through nervous systems. Basic to the subject is how single neurons 
transmit information. As in any modeling effort, we must ignore or summarize details and 
focus on what, we hope, are a few essential aspects. The ML model |20| has often been 
used as a good, qualitatively quite accurate, two-dimensional model of neuronal spiking. 
It is a conductance-based model like the HH model, introduced to explain the dynamics 
of the barnacle muscle fiber. The original ML model was three-dimensional, including a 
fast responding voltage-sensitive Ca^+ conductance, and a delayed voltage-dependent K+ 
conductance for recovery. To justify the two-dimensional version, one uses that the Ca^+ 
activation moves on a much faster time scale than the other variables, and can conveniently be 
treated as an instantaneous variable, by replacing it by its steady-state value given the other 
variables. 

The parameter values in our computations were chosen from Il22ll23l . and are given in Table 
[T]together with the interpretation of variables and parameters. The variable Vt represents the 
membrane potential of the neuron at time t, and Wt represents the normalized conductance of 
the K+ current. This is a variable between and 1, and could be interpreted as the probability 
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Table 1: Variables and parameter values used in the Morris-Lecar model 



V{t) [mV] Membrane voltage 

Vl^(i) [1] Normalized conductance 

t [ms] Time 







-1.2 mV 


Scaling parameter 


V2 




18 mV 


Scaling parameter 


^3 




2mV 


Scaling parameter 


Vi 




30 mV 


Scaling parameter 


5Ca 




4.4 ^S/cm^ 


Maximal conductance associated with Ca^+ current 


9K 




8 ^lSlcrc? 


Maximal conductance associated with K+ cuiTent 


9L 




2 ^S/cm^ 


Conductance associated with leak current 


Vc. 




120 mV 


Reversal potential for Ca^+ cuiTent 


Vk 




-84 mV 


Reversal potential for K+ current 


Vl 




-60 mV 


Reversal potential for leak current 


C 




20 ^F/cm^ 


Membrane capacitance 


4> 




0.04 1/ms 


Rate scaling parameter 


I 




90 ^A/cm^ 


Input current 



that a K+ ion channel is open at time t. The non-linear model equations are 

dVt - ^{-gcamoo(Vt){Vt-Vca)-gKWt{Vt-VK)-9L{Vt-VL)+I)dt, (1) 
dWt = {a{Vt){l-Wt)~ l3{Vt)Wt)dt, (2) 

with the auxiliary functions given by 

rriooiv) = i (l + tanh(^^^-^ ) ) , (3) 



2 V \ V-2 



"2 



aiv) ^ i0cosh(^)(l + tanh(!^)), (4) 
m = i^^coshf^Ul-tanhf^) ). (5) 



Equation ([T]l describing the dynamics of Vt contains four terms, coiTesponding to Ca^+ cuiTent, 
K+ current, a general leak current, and the input current /. The functions a( ) and /?(•) model 
the rates of opening and closing, respectively, of the K+ ion channels. The function moo( ) 
represents the equilibrium value of the normalized Ca^+ conductance for a given value of the 
membrane potential. 

In Fig. [TJ\ the phase-state of the model is plotted. The system has two stable attractors; 
a stable fixed point corresponding to quiescence of the neuron, and a stable limit cycle corre- 
sponding to repetitive firing. In between the two attractors is an unstable limit cycle, which 
splits the state space into two parts from either of which the deterministic process cannot 
escape, once trapped there. 

2.1. The stochastic Morris-Lecar model with channel noise 

It has long been known that the opening and closing of ion channels is an important part 
of neuron function. Channel activity is summarized, even in the comparatively detailed HH 
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Figure 1: Phase-state plots of the normalized conductance Wt against membrane voltage Vt- The full 
drawn magenta curve is a stable limit cycle, the dashed magenta curve is an unstable limit cycle, and the 
magenta point is a stable fixed point. Black curves are sample trajectories. Panel A; model without noise, 
(|T|-j5j. If the process is started between the stable and the unstable limit cycle, or outside the stable limit 
cycle, the solution is seen to spiral out, respectively in, towards the stable limit cycle, corresponding to 
repetitive firing of the neuron. If the process is started inside the unstable limit cycle, the solution spirals 
into the stable fixed point, corresponding to subthreshold fluctuations of the neuron. Note that three 
trajectories are plotted. Panel B: model with noise, fl}, ([3j-l[5j and ^, a* = 0.05. Only one trajectory 
is plotted, and the solution is seen to switch between periods of firing and quiescence. 

model, by potential dependent averages. However, it has become apparent that the stochastic 
nature of ion channels must be explicitly modeled if we are to capture essential features of 
neuron dynamics. Changes in the states of channels cannot be tracked explicitly because of 
their vast number. Hence, it is useful to model the role of ion channels as a stochastic process, 
Wt, the proportion of channels open at time t. We therefore add channel noise by changing 
the ordinary differential equation system ([T]) - (|5]), to a stochastic differential equation system, 
replacing the conductance equation (|2| by 

dWt = iaiVt)il-Wt)-l3iVt)Wt)dt + h{Vt,Wt)dBt, (6) 

where Bt is a standard Wiener process, and the function h{-) has to be chosen. 

The diffusion coefficient h{-) in ^ should be based on the drift coefficient which gives 
the rate of change of fraction of open ion channels due to random openings and closings. A 
natural choice of the function h{-), following the diffusion approximation of |16|, would be 
the square root of the sum of the two rates in the drift coefficient, times a factor 1 / \/N where 

is the number of ion channels involved. However, this choice has the problem that it is not 
zero when all the channels are closed, and the resulting (|6| would produce negative solutions 
with positive probability. To avoid this difficulty, for fixed Vt we let Wt be a Jacobi diffusion. 
In fact, in the class of Pearson diffusions |9 |, i.e. one-dimensional diffusions with linear drift, 
and with ) a polynomial of at most degree two, this is the only bounded diffusion. Living 
on (0, 1), it has the form 

dXt - -9 {Xt ~fi)dt + 7\/20Xt(l - Xt)dBt (7) 
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where 6 > Q and ji e (0, 1). It is named for the eigenfunctions of the generator, which are 
the Jacobi polynomials. It is ergodic provided that 7^ < min(/x, (1 — ^)), and its stationary 
distribution is the Beta distribution with shape parameters and (1 — /i)/7^. It has mean /i 
and variance 7^/Lt(l — /i)/(l + 7^). In our case, because the diffusion coefficient in (|7]) should 
be of the same order as the one given by the Kurtz approximation IIT6I . 7 is proportional to 

By equating the drift terms in (|6]l and (|7|, we have 6 = a{Vt) + l5{Vt) and /i = a{Vt)/ 
{a{Vt)+P{Vt)). So for fixed Vt, with h-'iVt, Wt) = -f^2{a(Vt)+p{Vt))Wt[l-Wt), where 7^ 
is constrained by 7^ (a(Vt)+/3(Vt)) < min(Q;(Vt), /3(Vt)), also (|6| will stay bounded in (0, 1). 
Since a{Vt) and are strictly positive, we can put 7^ = {(j*Yaiyt)j3{Vt)/{a{Vt) + 

/3(Vf))^, with cr* e (0, 1], and specify the conductance equation (|6]l as 

dWt = {a{Vt){l~Wt)-p{Vt)Wt)dt + G*J2^^^^^^^Wt{l-Wt)dBt.{%) 

In the next Section we compute the equilibrium point (Veq, Weq) of the system ([T]i-(|5]l for 
the chosen parameters. By equating the diffusion coefficient as it would occur in the diffusion 
approximation of |,16| with the one in ([H) at (Veq, Weq) we will obtain a* in terms of 1/\/N, 
where N is the number of channels involved. 

It can be shown by a coupling argument that also for varying Vt will Wt given by ^ stay 
bounded in (0, 1), since Vt is bounded once it is started inside some interval [Tj. 

In Fig. |2j the model defined by ([T]l, ([3]l-(|5]) and ([HJ is simulated for different values of a*, 
where these can be thought of as corresponding to different total numbers of ion channels. 



3. The linear approximation of the stochastic Morris-Lecar during quiescence 

To identify the process of subthreshold oscillations, i.e. the dynamics close to the stable 
fixed point between firings, we analyze the linearized system around this point. Consider the 
system 

dVt = f{Vt,Wt)dt, 
dWt = g{Vt,Wt)dt + h{Vt,Wt)dBt, 

where the functions /(•), g{-) and h{-) are given by ([T]i, (|3]l-(|5]l and ([H). 

For the chosen parameter values given in Table [T] the deterministic system, obtained for 
h{-) = 0, has a unique locally stable equilibrium point (T4q, Weq) given by 



= a(V.,)Vff(14,) '2l.^ + '-°"v V, 

and V^q is the solution to the equation /(V^q, Weq(V;:q)) — 0, which cannot be solved analyti- 
cally, but can be found numerically. The input current value / = 90/iA/cm^ is a typical value 
well inside the range of / where the deterministic dynamics has a stable limit point inside an 
unstable limit cycle as shown in Fig. The equilibrium point for / = 90/iA/cm^ is 

(V;q,Weq) - (-26.6 mV , 0.129). 

In terms of the centered variables 

Xi^^ =Vt- 14q , ^ ^Wt- Weq 
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Figure 2: Time series plots (black curves) of the stochastic Morris-Lecar model for different noise 
levels started inside the unstable limit cycle, but not at the fixed point. Upper left: a* = 0.02, upper 
right: a* — 0.03, lower left: a* = 0.05, lower right: a* — 0.1. Note different scales, in the upper left 
panel there is no firing. The magenta curves are the deterministic model, a* = 0. 



the system becomes 



(1) 



(2) 



dX, 



We write X* 



nx['\x[^^)dt, 



(2) 



+ l^eq)) dt + h + Feq), (^f ^ + Weq)) ) 



(9) 



h*{Xl^\X, 



(1) 



)dB^\ 



(10) 



Note that b'"^^ does not enter the dynamics, but is introduced to ease the matrix notation, as 
will be clear in the following. When the noise is small and the process Xt is started near the 
equilibrium point, x — (0, 0), we expect the dynamics to concentrate around the equilibrium 
point. A local approximation is obtained by linearizing (|9ll-(fT0| around (0, 0). The diffusion 



xl^^)^ and Bt ~ {b[^^ B\^'')'^, where denotes transposition. 



3(2)^T 
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term is approximated by setting x[^^ = x'f'^ = in the diffusion coefficients. The linearized 
system is 

dXt = mXtdt + GdBu (11) 



where 
M 



f mil 


■mi2\ 




(df 






1 ^^i 
I dal 


dx2 1 


\m2i 


TO22/ 




9X2 / 








\dxi 



(xi,a;2) = (0,0) 



0.0258 -22.96l\ 
0.000335 -0.0446 J 



using the parameter values in Table [T] and 




/O 

VO ^V2(a(Kq) +/3(Kq))(l - Weq)T¥ec ^ " ' ^^^^ 

where a — 0.034cr*. By evaluating the diffusion approximation of lfT6l at (Kq, Weq) and 
equating to the above we obtain a* = l/-\/Weq(l — Weq)N 3/VN. In the Appendix the 
matrix M is detailed. Solutions of ( [TT] i with G = are given in terms of the eigenvalues of 
M which are complex conjugates and given by 

-X±uji = -0.0094 ± 0.0803i 

where A — — tr(M)/2, = — det(M)| and i = \/— T. Thus, near the equilibrium point 
the solution of ( fTT) , with ct — 0, is 

^ ch'^f)e-^\ (13) 

where C contains the initial conditions 

^xq (mi2?/o + {mil + X)xa)/uj 



C — 

\yo {m2iXo - {mil + X)ya)/uj 

In Fig. |3]the solution of the deterministic model, ([T])-(|5]l with ct = 0, is compared to the linear 
approximation ( [T3] l. 



4. Identification of the stochastic process of quiescence 

In this Section we identify the stochastic process defined by the linearized system ( fTT) in the 
limit of small A, i.e. under the condition X u. The deterministic system ( [T3| ) has decaying 
oscillations, whereas for the stochastic system ( fTT) , the noise will prevent the decay of the 
oscillations. Can we describe the resulting process specifically? The answer is that, after a 
linear change of variables, this process can be approximated in distribution by a fixed matrix 
times a deterministic circular motion modulated by an OU process. 

We follow the development in [2l, where a first step is to transform the matrix M into a form 
which reveals the slow decay towards the equilibrium point and the fast oscillatory structure 
of the deterministic dynamics. Let Q be a 2 x 2 matrix such that 

Q-'MQ . . A, 
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Figure 3: The solution of the 
deterministic model ([T]|-(|5]l with 
a = (black full drawn curves) 
is compared to the linear approxi- 
mation ( [T3| l (cyan dashed curves). 
Upper panel: normalized conduc- 
tance Wt (dimensionless). Lower 
panel: membrane potential Vt 
(mV). Time is measured in ms. 
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A possible choice for Q is 



Leti:* = Q^^Xt, then 



Q 



—u! mil + A 

17121 



dXt = AXtdt + CdBt 



(14) 



where C = Q ^G. A further change of variables moves the rotation to form part of the 
diffusion coefficient of the linear stochastic system. We define 



where 



Rs = 



RujfXt 



cos s — sm s 



Vsins cos s 

is the counterclockwise rotation of angle s. Then by Ito's formula 

dXt = -\Xtdt + R^tCdBf 
The infinitesimal covariance matrix in ([14]) is 

B = CC^ ^ Q-iGG^(Q-i)^ = ( ^7" + ^!' -(-11 + A) 

Now define 



^tr(B) = \{Bii+B22) 



a'^mi2 
2uj^m2i ' 



(15) 



(16) 
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where we have used that (mn + A)^ -\- uj"^ — — 1711271121- Finally, we rescale Xt so that we 
can compare with a standardized two-dimensional OU process. Let 



Relation ([TSj becomes 



Ut — Xf/x- 



dUt = -Utdt+^R^t/),CdBt (17) 



where Bt — V^Bf /x is another standard two-dimensional Brownian motion. The following 
Theorem from |2j allows us to approximate the process Ut given by ([TtJ, by a two-dimensional 
OU process with independent coordinates. 

Theorem 4.1. For each fixed t* > Q and x e the distribution of {Ut : < t < t*} 
given by with Uq = x converges as X/cu — > to the distribution of the standardized 
two-dimensional OU process {St : < t < t*} generated by 

dSt = -Stdt + dBt 

with Sq = X. 

Here St follows a normal distribution, St N (^qg"*, ^(1 — e~^*)I), where I is the 2 x 2 
identity matrix. The proof of this Theorem uses a martingale problem convergence argument 
and involves the notion of stochastic averaging, where fast oscillations integrate out revealing 
the remaining structure determined by slower oscillations. Another result of this type obtained 
by a different method, called multiscale analysis, is in lUTll . 

Thus, the process Ut is approximated by 5t if A ^ uj. In our case A is one order of 
magnitude smaller than lu. 

Putting together the transformations and the final approximation we have, in the sense of 
stochastic process distributions, 

~ ~ T T 

Xt = QXt = QR^uitXt = QR-uit^^U\t ~ QR-LJt—^Sxt 



Sxf (18) 



T /— mil + A\ f cosLut sinutt 
\ 77121 / I — sinwt cosut 

Let us denote by Xt the stochastic process on the right hand side of ([TSj, i.e. 

X^ = TQR_^tSxtl^- (19) 

To get a sense of how closely the process Xt approximates the dynamics of the ML process in 
a neighborhood of (T4q, Weq) we compare their power spectral densities, as well as that of the 
solution of the linearized system ( [TT| . The spectral density of X^ and that ofXt satisfying ([TTJ 
can be calculated explicitly using the power spectrum formula of 1 10| for linear diffusions of 
the form ( fTTj l. In fact Xf is such a diffusion: the effect of the stochastic averaging can be seen 
as replacing C from ( [l4j l by a multiple of the identity in the system ( [l4j i, so the approximation 
to X satisfies dXf ~ AXtdt + rdBt, where r is given by ( fTS) !. If we transform this equation 
by Xt = QXt, we see that X" satisfies 

dX^ = MX^dt + rQdBt. (20) 
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The spectral density of the first coordinate of X°- is 

^ (/^+det(M)) 

27r ((/2 - det(M))2 + (/tr(M))2) 2cj2 

whereas the spectral density of the first coordinate of the linearized system, ( [TT] i, is 



27r (/2 - det(M))2 + (/tr(M))2 



In Fig. |4]the theoretical spectral densities for the two approximations are plotted, together 
with the estimated spectral density of the quiescent process from simulations of the stochastic 
ML model ([T]l,(|3]l-(|5]l and (|8]l. The spectral density is estimated by averaging over at least 20 
estimates from paths started at of at least 450 ms of subthreshold fluctuations, and scaled to 
have the same maximum as the theoretical spectral density from ( |20] l. The averaging is done to 
reduce the large variance connected with spectral density estimation, avoiding any smoothing. 
Thus, the estimator is approximately unbiased, see also 1 8 1 where this approach is treated. The 
estimation is done for <j* — 0.03, 0.05 and 0.1. For higher noise, the lengths of subthreshold 
fluctuations between spikes are too short to reliably estimate the spectral density. Moreover, 
a* =0.1 corresponds to a number of ion channels N k, 900, which can be considered a 
minimum acceptable number for the diffusion approximation to be relevant. The value <j* = 
0.03 corresponds to iV sa 10, 000. Remember that a = 0.034cr*, see ( fT2] i. 

The approximations are only acceptable for small noise, which is expected, since larger 
noise brings the process to areas further away from the fixed point, where non-linearities 
become increasingly important. 




0.0 (0 0.1 0.2 0.0 (0 0.1 0.2 0.0 (0 0.1 0.2 

frequency [1/ms] frequency [1/ms] frequency [1/ms] 



Figure 4: Spectral density estimated from simulations between spikes of model ([TJ, j3j-([5}, ([8} (black 
solid line), theoretical spectral density of model \2Q) (cyan dashed line), and theoretical spectral density 
of model \\\) (magenta dotted line). From left to right: cr* = 0.03, 0.05 and 0.1. 



5. Reconstructing the stochastic ML firing mechanism 

In this Section we construct a firing mechanism matching that of the stochastic ML neuron. 
In Section|6]we will define a new LIF-type process by combining this firing mechanism with 
the radial OU process. This new model will, for small cr, have an ISI distribution similar to 
that of the ML. 
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Firing in model ([TJ, ([3]l-(|5]l and ([8]) occurs when the stochastic dynamics shifts from a path 
circulating the stable equilibrium, modulated by an OU, to a noisy circuiting of the stable 
limit cycle. This shift happens, roughly, when the orbit passes from the inside to the outside 
of the unstable limit cycle. When the orbit comes close to the unstable limit cycle, it will 
follow this limit cycle for a short time, and then escape either to the inside, i.e. continue 
its subthreshold oscillations, or to the outside and a spike will occur. This understanding is 
not accurate enough to be implemented as a firing scheme for the radial OU process ( p7j l, 
as we discuss further in Section |7] Hence, we embed the process X" defined by ([T9| in the 
stochastic ML model by constructing a firing mechanism mimicking that of the ML itself. 
It is clear that in the ML model, starting inside the unstable limit cycle, a spike will occur 
with increasing probability, the further away the process is from the fixed point. In order to 
construct a firing mechanism matching that of ML, we will estimate, from simulations, the 
conditional probability that the ML fires, given that the trajectory of the ML crosses the line 
L = {{v^w) : V = Veq,w < Weq}. We computed estimates from simulated data using 
crossings of the line L as follows. 

For a given value of a* and distance I from the fixed point, a short trajectory starting in 
(Veq, Weq — 1) was simulated from model ([TJ, (|3]l-(|5]l and ([Sj, and it was registered whether 
firing occurred in the first cycle of the stochastic path around (Veq, W^eq)- Firing was defined 
by the path crossing the line v = 0, which is well above the largest level inside the unstable 
limit cycle, see Fig. [T]3. This was repeated 1000 times, and estimates of the conditional 
probability of spiking, p{l,(T*), were computed as the frequency of the trajectories where 
firing occurred. The procedure was repeated for / — U = i5,i — 1, . . . , 25, where 6 is the 
distance to the stable limit cycle divided by 20. In this way a grid of possible I values was 
covered, starting from ^ = at the fixed point, where the probability of firing is close to zero, 
to a point on L below the stable limit cycle, where the probability of firing is close to one. The 
estimation was, furthermore, repeated for a* = 0.01 to 0.08 in steps of 0.01. 

For each fixed a* , the estimates of the conditional probability appear to depend in a sig- 
moidal way on the distance from the fixed point. We assumed the conditional firing probability 
to be of the form 

1 + exp((a - 1)1 p) 

The parameters a and were estimated using non-linear regression of the 25 estimates of 
p{li;a*) on I. In Fig. |5|\ these parametric estimates are plotted, as well as the individual 
nonparametric estimates p for a* = 0.02, 0.05 and 0.08. We see that the family of estimates, 
p, fits the hypothetic curve quite well for each value of a* . Regression estimates are reported 
in Table |2] Note that a is the distance along L from Weq at which the conditional probability 
of firing equals one half. For all values of cr*, the estimate of a is close to the distance along L 
between T^^eq and the unstable limit cycle, which equals 0.0172. In other words, the probability 



cr* 


0.01 


0.02 


0.03 


0.04 


0.05 


0.06 


0.07 


0.08 


a 


0.0174 


0.0174 


0.0169 


0.0168 


0.0171 


0.0169 


0.0167 


0.0168 




0.0006 


0.0013 


0.0020 


0.0028 


0.0033 


0.0039 


0.0047 


0.0054 




7.1022 


3.5426 


2.3012 


1.7156 


1.3922 


1.1474 


0.9739 


0.8549 




0.2590 


0.2624 


0.2759 


0.2831 


0.2718 


0.2674 


0.2738 


0.2764 



Table 2: Estimates of regression parameters for p(-) in the original space (first two rows), and in the 
transformed coordinates (last two rows). 
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0.005 0.015 0.025 1 2 3 4 

distance from fixed point distance in transformed space 



Figure 5: Conditional probability of spiking when crossing the line L — {{Vjw) : v — l^eq,'!^' < 
Weq} for different values of a*. (A) Original space. The circles, plus's and stars are individual 
nonparametric estimates obtained using a* = 0.02, 0.05 and 0.08, respectively, with the fitted curves 
on top given by |2TJ. The dashed line indicates where the unstable limit cycle crosses L, the full 
drawn line where the stable limit cycle crosses L. (B) The fitted curves in the transformed space for 
a* = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07 and 0.08 (right to left), as a function of the distance from the 
fixed point in the transformed coordinates. The crosses and boxed crosses indicate the crossing of the 
unstable and stable limit cycles of L, respectively, which depend on a = 0.034cr*. 

of firing, if the path starts at the intersection of L with the unstable limit cycle, is about 1/2. 
The parameter /3 indicates the width of a band around a where the conditional probability 
essentially changes. For instance, if I E a ± /3 then p{l) e (0.27, 0.73), if I £ a zL 2f3 then 
p{l) € (0.12,0.88). As expected, the estimate of /3 increases with increasing a*, and for 
small noise the conditional probability approaches a step function since the process is mostly 
dominated by the drift. A step function would correspond to the firing being represented by 
a first-passage time of a fixed threshold. Note though that (3 is approximately proportional to 
a*, and thus, as we said earlier and will see in the following, a fixed threshold at the crossing 
of the unstable limit cycle does not reproduce the desired spiking characteristics. 

In order to simplify the construction in Section |6] of a LIF model which, together with a 
firing rule, behaves like the stochastic ML, we will change coordinates as follows. Observe 
that ^T9\ can be written 

^Q-^Xf = f^"^"^^ (22) 
r V— smojt cos u}t I 



so for fixed t, -\/AQ ^X^ /t is the clockwise rotation by angle ut of the orthogonal pair 
{S[\\S''^^^). We define a tr 
normalizing as in (|22]l. Let 



(S'^j'', S'^t''). We define a transformation of the space (w, w) by centering at (Veq, Weq) and 



~ , Q"M w (23) 

Wj T \W- Weq/ 



be the coordinates of the transformed space. In the new coordinates our process is simplified 
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to a rotation modulated by a standard two-dimensional OU process with independent compo- 
nents. 

The transformation depends on cr = 0.034cr*, namely, the transformed unstable limit cycle 
becomes smaller with increasing noise, through the value of r given in ( [T6] l. This is exactly 
what is causing a higher firing probability for larger a* . The line L will in the transformed 
space be 

T \l I TO2ir \ 1 / 



for ^ > 0. A distance I will thus transform to a distance r — {\/2\/ij)l, and the conditional 
probability of firing pT| ) transforms to 

1 -I- exp((a* — r)/p*j 



where a* = aV2X/a and /3* = l3V2X/a. The fitted curves of (|24| for a* = 0.02 - 0.08, 
as a function of the distance from the fixed point in the transformed coordinates are given in 
Fig. |5|3, with indication of the crossings of the unstable and stable limit cycles, respectively, 
which now depend on a. Note that in the transformed space, the width of the band where the 
conditional probability is essentially different from or 1 is nearly constant, see Table[2] From 
here on we use the coordinates defined by (|23]l. 



6. Construction of a leaky-integrate-and-fire model with ML firing statistics 

The simpler stochastic LIF models sacrifice realism for mathematical tractability |l4l[TT|. In 
these models, a neuron is characterized by a single stochastic differential equation describing 
the evolution of neuronal membrane potential depending on time, 

dXt = iiiXt)dt + a{Xt)dBt, Xo=xo, (25) 

where Xt corresponds to Vt in the ML model, together with a threshold firing rule, 

T = M{t >0:Xt>S}. (26) 

In this Section we define a LIF model which does not make this compromise, using the result 
of Section|4]and the firing mechanism defined in Section|5] 

The distance of the approximate process y/XQ^^X^ /t of ( p2] i from the point (0, 0) at time 
t is given by the modulus of the two-dimensional standardized OU process S\t- The modulus 
of Sxt at time t is given by the process 



Rxt = 



which is a standard radial OU process with two degrees of freedom. It has state space (0, oo), 
and solves the stochastic differential equation 

dRxt = (2^^^^^*)'^^ + ^^^*' ^27) 

see e.g. HJ. We define a new LIF process by p7| ), and firing mechanism derived from ( |24| i. 
After each firing, we will reset the time to and assume the process reset to 0, i.e. Rq — 0, 
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corresponding to 5*0 = (0,0) and (Vo,VFo) = (Veq,Weq)- By Ito's formula, the process 
y„ = satisfies the stochastic differential equation 

dYu = 2{l-Yu)du + 2y/Y^dWu, (28) 

and is thus a square-root process, see e.g. [5|, also called a Feller or a Cox-Ingersoll-Ross 
process. This process is ergodic, and its stationary distribution is the exponential distribution 
with mean one. It follows that the stationary distribution of _R„ has density /(r) = 2re^^ on 
(0, oo), i.e. it follows a Rayleigh distribution. The transition density of y„ starting at j/o at time 
0, is a non-central -distribution with two degrees of freedom and non-centrality parameter 
S{u, uo) — 2j/oe^^"/(l ^ e^^"). Then 2y„/(l — e^^") follows the standard non-central x^- 
distribution F^2{2y/{1 — e^^"), 2,6{u,yo)). It is particularly simple because of the integer 
degrees of freedom. Transforming to the radial OU we obtain the transition density of i?„ 
starting at s at time 

2r f + s^e~^" 1 f \ 

U-^^) = T3^-P|-^— ^l^ol^-j^j, (29) 

where Iq {x) = ^ ^°^^dO is the modified Bessel function of the first kind of index 0. 

Writing the two-dimensional process Su in polar coordinates, i?„ and 9u, where 6u is the 
angle at time u to the positive part of the first coordinate, we find that the modulus and the 
angle are independent, and that 6^ is uniformly distributed on (0, 2tt). This can e.g. be seen 
from the fact that S'i^'' and S'i'^'' are independent normal with mean and equal variances. 
Thus, for fixed u, si^^ /Su'^ is standard Cauchy distributed and 0„ = arctan(S'i'^''/S'i^'') is 
C/(0,27r). 

Let T denote the firing time random variable. We want to compute the density of the 
distribution of T, and for this we find it convenient to express this density in terms of the 
conditional hazard rate, 

a(t,r) = liiR ^P{t<T <t + At\T>t,Rxt=r). 

At^Q At 

This function is the density of the conditional probability, given the position on i is r at time 
t, of a spike occuring in the next small time interval, given that it has not yet occurred. 
From standard results from survival analysis, see e.g. 1 1 1, we obtain 



PiT >t\Rxs,0<s<t) = exp^-^a 



The unconditional distribution is then given by 



{Rxs)ds 



a{Rxs)ds (30) 



P{T >t) ^ E ^exp 
where E{-) denotes expectation with respect to the distribution of R. The density is thus 

g{t) = j^P{T <t) = E(^aiRxt)c^p(^- J^a{Rxs)ds^y (31) 
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The firing is defined to be initiated from L, and on average the process crosses L every 2'k/uj — 
78.2 time units. Using ( p4| i, the estimated conditional probabiUty of firing given the position 
on L is r, which by definition does not depend on t, we estimate the hazard rate as 

a{t,r) = air) = — rrTTT- (32) 

^ ' ' ^ ' 2tt 1 + exp((a* - r)//3*) 

Note that it is bounded. This is not realistic, since a very large value of r should cause 
immediate firing. 

In ETi a firing rule with unbounded hazard rate was proposed, and in [151 it was shown to 
fit well to experimental data. Therefore, we will also see how our model performs if we use in 
the firing mechanism a hazard rate of the form 

a{t,r) — a[r) = exp((r — a)//3) (33) 

for a, /? > 0. Like before, a plays the role of a threshold, and j3 gives the width of the threshold 
region. When f3 ^ Q, the firing rule converges to a fixed threshold crossing. To estimate a and 
/3 in ( (33] l, we simulated 1000 spike times from the ML. The cumulative hazard A{t) — a{t) 
was then estimated from the simulated spike times by the standard empirical Nelson-Aalen 
estimator. The theoretical cumulative hazard using (|27]i and ((33]l can be calculated as 



A{t) = ^ a{Rxs)ds^ = exp (^-^ j^E (^exp ds 

= V^exp(^-^^J^ (^.g(5)expQg(,s)2^ $(g(,s)) + l^d5 (34) 

where we have used the density f\s{r, 0) given in (|29]). Here, g{s) = W^^^P^/P, and 
$(•) is the standard normal cumulative distribution function. Then, a and /3 were estimated by 
the least square distance between ( (34] i and the estimated cumulative hazard from the simulated 
spike times. For cr* — 0.05 the estimates were a = 6.31 and (3 = 0.76. 
The final model is 

dRu - {^~Ru]dt + dWu-Ru~lj^{Ru-.du), (35) 



2Ru 

where iJ,{Ru~ , du) is a Poisson measure with intensity Q;(i?„_ ), and denotes the left limit 
of Ru- Here, a(-) is either given by ( (32] i or ( (33] l. The jump size is —R^-, thus giving the reset 
to at spike times. 

A reasonable alternative to the soft threshold firing mechanism used here would be to use the 
firing rule defined by a threshold as in the classical LIF models, equation (|26|. A natural choice 
of threshold would be where the LIF process reaches a level corresponding to the unstable limit 
cycle. In fact, according to our estimates in Fig. [5] and Table |2] the firing probability of the 
ML at this threshold is around 1/2. However, the ISI distribution estimated from simulations 
using a hard threshold at the unstable limit cycle is shifted towards larger times, relative to 
the ML ISI distribution. This happens because the process might cycle many times inside the 
unstable limit cycle, so even if the probability of spiking in a single cycle is small, the total 
probability is not negligible. This is lost when only a hard threshold is considered. Instead we 
chose the threshold value such that the mean of the ML ISI distribution and the mean of the 
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LIF ISI distribution were the same. In [121, the mean of T from (|26]l with Xt — Rt started at 
i?o = is given using a hypergeometric function, 

E{T) = ^ 2^2 (1,1; 2, 2; ^2)^ (36) 

The average of the 1000 ML firing times for a* — 0.05 was 447. Equating with ( (36] l gives 
a value S = 2.97 for the hard threshold. Note that this is much smaller than the estimated a 
from 

7. Comparison of firing statistics 

One of the major issues in computational neuroscience is to determine the ISI distribution. 
We therefore simulated the ML model given by ([TJ, (|3]l-(|5]l and ([8]) until spiking, and thereafter 
reset to the fixed point. This was done 1000 times, and the time of the firing was recorded. 
The ISI distribution from our approximate model is given by the density ( (3T| , or equivalently, 
from the survival function ([30|. Due to the law of large numbers and since we know the exact 
distribution of i?„, for fixed t we can numerically determine ( (3T] l up to any desired precision 
by choosing n and M large enough through the expression 

m— 1 \ i—1 I 

Here (i?o™^ ■ • • ' -^iAt/n' ■ • ■ ' ^aT^) ^ realizations of Ri\t/n^ i — 0, 1, . . . , n, and the 
integral has been approximated by the trapezoidal rule. The hazard rate is either given by ( (32] i 
or(|33]). 

The results are illustrated in Figure |6]for a* = 0.05, using M — 1000. The estimated ISI 
distributions from our approximate model with both firing mechanisms compare well with the 
estimated ISI distribution of ML reset to after firings. On the contrary, the hard threshold 
does not reproduce the ISI distribution well, e.g. the right tail is too heavy. This is because the 
probability of firing during low subthreshold activity is set to 0, whereas we have seen it is not. 

8. Discussion 

A stochastic LIF model constructed with a radial OU process and firing mechanism of either 
logistic or exponential type has been shown to mimic the ISI statistics of a ML neuron model. 
It captures subthreshold dynamics, not of the membrane potential alone, but of a combination 
of the membrane potential and ion channels. This construction will allow us to answer several 
questions about ML models, which have been accessible only for LIF models, even though the 
latter have less biological motivation. 

An example of such a question would be: Using ISI experimental data, the noise standard 
deviation a can be estimated 1 18 1. In principle, this should also be possible from our new LIF 
model, even though we use a soft threshold. This will give an estimate of N, the number of 
ion channels involved, through the relation (cr*)^ ~ 

A question we have not explored is: what is the best way to restart our new LIF model? In 
our simulations we restarted both our LIF and the ML at the fixed point of the ML. However, 
an uninterrupted stochastic ML produces continuous paths as in Fig. [TJi. After firing, which 
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Figure 6: Distribution of firing 
times for a* = 0.05. The his- 
togram is based on 1000 simulated 
firing times from the ML model, 
the vertical dotted line is the aver- 
age. Curves are estimates of the 
probability density, equation ( |37| i. 
Black curve is estimated using 
32 1, gray curve is estimated using 
33ll, dashed curve is estimated 



using a fixed threshold, d26| 



means traversing the large stable limit cycle, possibly several times, they reenter a neighbor- 
hood of the fixed point from its edge. A further refinement of our LIF model will be obtained 
by introducing a reentry mechanism, which mimics this aspect of the ML. 



Appendix A. Linearization matrix 

The expression for M in ( [TT| i is 



M 



' mn -gfeWeq(T/eq - Vk)/C\ 

^214qW^eq/3 (V,^) /V4 -a (Kq) / 



^ 2ffCa(T4q-^Ca)a(Kq)/?(Kq) _^ . . ^ „^ ^ 
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